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Abstract. How, in principle, could one solve the atomic 
structure of a quasicrystal, modeled as a random tiling dec- 
orated by atoms, and what techniques are available to do 
it? One path is to solve the phase problem first, obtaining 
the density in a higher dimensional space which yields the 
averaged scattering density in 3-dimensional space by the 
usual construction of an incommensurate cut. A novel direct 
method for this is summarized and applied to an i(AlPdMn) 
data set. This averaged density falls short of a true structure 
(ietermination (which would reveal the typical unaveraged 
atomic patterns.) We discuss the problematic validity of in- 
ferring an ideal structure by simply factoring out a "perp- 
space" Debye- Waller factor, and we test this using simula- 
tions of rhombohedral tilings. A second, "unified" path is to 
relate the measured and modeled intensities directly, by ad- 
justing parameters in a simulation to optimize the fit. This 
approach is well suited for unifying structural information 
from diffraction and from minimizing total energies derived 
ultimately from ab-initio calculations. Finally, we discuss 
the special pitfalls of fitting random- tiling decagonal phases. 



1 Introduction 

The "random tiling" model of quasicrystals has specific 
and sometimes radical implications for the procedures that 
should be used to determine the atomic structures. This pa- 
per collects several kinds of answer to the question, "How 
should one modify the standard crystallographic procedures 
in order to solve the atomic structure of a random tiling 
quasicrystal?" (Notice that even the definition of "solve" is 
arguable in the case of an intrinsically random structure.) 
We intend this paper to encourage crystallographers to try 
random-tiling fits to diffraction data, by showing that a cou- 
ple of approaches have already been thought out in some de- 
tail. It also cautions against some common misunderstand- 
ings about this point, and collects three new (and imperfectly 
digested) numerical tests of our ideas. The discussion goes 
well beyond its antecedent in Sec. 7 of Ref. 
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1.1 Perp space and diffraction 

The "perp" coordinate is a central notion in quasicrys- 
tallography. Before proceeding, we had better remind the 
reader that it has two definitions, which coincide for the case 
of an ideal tiling. 

Definition 1, valid for the usual higher-dimensional 
"cut" description [Q] of any quasiperiodic struc- 
ture (e.g. an ideal tiling). The (periodic) higher- 
dimensional density consists of atomic surfaces ex- 
tending in the perp direction; then is the perp- 
space displacement from the intersection of the cut 
plane and the atomic surface, to a vertex of the higher- 
dimensional lattice. 

Definition 2, valid for any arbitrary (e.g. random) 
tiling. Any tiling can be "lifted" so that each tile ver- 
tex maps to a vertex of a periodic lattice in a higher 
dimension [[l|, ^; then of each tile vertex is the 
difference between the perp component of its lattice 
vector and the mean value ("center-of-mass" in perp 
space). 

Consider the probabihty distribution of the Definition 
2 h-^ values. It is easy to check that, as long as the physical 
(hyper)plane has an irrational orientation, and in the limit of 
a large physical-space volume V, there are ni^yd^h.^ possi- 
ble perp-space points in a perp volume element d^h.^, where 
Uq^^ is the volume of the unit cell in the -dimensional lat- 
tice. The mean number of such points actually present in the 
tihng is defined to be p{h.^)naV d^h.^ , so p{h.-^) e [0, 1] is 
a probability. In random tilings it is, in practice, a smooth 
function. (In the special case of an ideal tiling, pih.^) is 1 
within the "acceptance domain" where sites are sure to be 
occupied, and zero outside it.) Note h-'", according to Defi- 
nition 1, or its mean h-'- = J d^h.^h.^ p{h.^) according to 
Definition 2, plays the role for perp displacements that a 
crystal's center-of-mass plays for ordinary displacements. 

Now, we may take an ensemble average of a ran- 
dom tiling ensemble, just constraining the value to be 
zero (without loss of generality). This gives, in real space, a 
nonzero density p{y) on a discrete (but dense) set of points 
r. It is a fact, natural though not rigorously proven, that p{y) 
depends only on the perp-coordinate of the point r. If 
so, it is easy to see p(r) — pih.^) as defined in the last para- 
graph. In other words, starting with Definition 2, we found 
the averaged density p(r) is given by the cut construction of 
Definition I, if the atomic surfaces are weighted by p{h.^. 
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This shows that p(r) is quasiperiodic, since that term is, 
practically speaking, defined by the existence of a cut con- 
struction. [] 

These facts about pih.^) easily generalize to atomic 
decorations in which atoms sit on some of the tile vertices. 
Thus wherever the physical 3-plane cuts an atomic surface, 
p(h^) is the probability that an atom is present. 

The Bragg component of the diffraction amplitude Fg 
is defined as the part which scales as V with the volume. 
It is, in fact, simply the Fourier transform of the averaged 
density: 



(1) 



Here is the perp-space partner of the reciprocal lat- 
tice wavevector g. In cases other than icosahedral, "3" in 
eq. is replaced by "d±", the dimensionality of perp space. 
(However we shall assume the physical dimension is always 
3, since this paper considers determinations from Bragg 
peaks, but random tilings have no Bragg peaks in dimen- 
sion 1 or 2.) For a tiling decorated by real atoms, (|l]) can be 
suitably generalized by including form factors and allowing 
for different classes of site (see eq. (^7|), below). 

1.2 Random tilings and the cut description 

Bragg diffraction amplitudes are always represented by a 
quasiperiodic density in higher-dimensional space, which 
when cut by a 3-plane at the correct incommensurate ori- 
entation gives the diffracting density in physical space. [||] 
The crystallographic refinements done to date have presup- 
posed that the structure, ideally, is quasiperiodic. But an (at 
least) equally plausible scenario is that the thermal equilib- 
rium structure is an intrinsically random tiling and so that its 
long-range quasiperiodic order is stabilized by the tiling's 
configurational entropy. |l]] Notions concerning the free 
energy and elastic constants are central to the physics [^, |l]] 
but will not be repeated here since they are somewhat tan- 
gential to the crystallography. Indeed, large parts of this 
paper could be applied equally well to the cases of (i) the 
"icosahedral glass" [^, e.g. i(AlCuLi), or i(TiNiZr), which 
is not thermodynamically stable; or (ii) the "weak matching- 
rule" quasicrystal which is a slight generalization of the 
ideal, energetically stabilized quasicrystal that nevertheless 
has a nonzero density of intrinsically random sites. 

So, how does the cut description get modified in the 
case of a random tiling? As just noted in Sec. 1.1, the 



physical-space density giving rise to the Bragg amplitudes is 
the ensemble average of the scattering density. This average 
is perfectly quasiperiodic, and so even for a random tiling it 
is represented as a cut by a 3-plane through a function in a 
higher-dimensional space. 

However, the physical-space density will have a great 
many fractionally occupied atoms; even if there is no sub- 
stitutional disorder, since a given volume of space has prob- 
abilities to be divided into tiles in different ways and dif- 
ferent tiles have different decorations. A second effect is 



In the case of the two-dimensional equilibrium random tiling or the 
icosahedral glass, the above arguments break down because p(h ) — 
0, i.e. the distribution gets broader and broader with increasing V. 



the "physical-space displacements" the deviations of an 
atom from ideal tiling vertex positions in response to forces 
from neighboring atoms. In ideal structures, an atom's lo- 
cal environment is a deterministic function of . Then the 
atom's equilibrium position is Tq + u(h^), also a function 
of h^, and this is the parametrization of the atomic surface's 
shape. In the random case, however, atoms corresponding 
to the same can have different local environments and 
hence different parallel-space displacements. This manifests 
itself as split positions in the physical-space cut, or equiva- 
lently the atomic surface becomes an discrete family of sur- 
faces, displaced from each other by small offsets in parallel 
space, and each having its own distribution p{h^)- 

Crystallographers would sometimes argue that a 
Fourier map reconstructing the correct ensemble-average 
density is "the structure" and is the proper goal when we 
fit the diffraction. In other words, they would say the task is 
done as soon as the phase problem is solved. We disagree 
strongly with this viewpoint. The aim of structural studies is 
to uncover the actual structure, which means understanding 
what a typical realization is like and not just the average over 
realizations. If the averaged structure contains two nearby, 
half-occupied sites, does the actual structure contain exactly 
one atom which can occupy either site at random? Or is it 
that 50% of the time both sites are occupied, and at other 
times some other site is occupied? Spurious but rare sites 
have negligible effect on the i?-f actor, but major effects on 
the total energy and the electronic states computed from the 
structure. Hence, the goal should be a proper description of 
the entire ensemble of configurations. 

1.3 Outline of the paper 

There are two alternate approaches for solving a random 
tiling structure: 

1. "Standard" approach. Break up the task into two 
stages. The first stage is to solve the phase problem with- 
out a complete model of the real-space structure. The (aver- 
aged) scattering density p(r) can then be constructed as an 
incommensurate cut in a high-dimensional space. The sec- 
ond stage is to relate this density to an atomic model. 

2. "Unified" approach. Assume a particular kind of 
random-tiling model, with variable parameters (such as the 
positions and species of decorating atoms on the tiles, or co- 
efficients in the "tile Hamiltonian" that governs the statistics 
of patterns in the random tiling). Simulate the tiling, calcu- 
late the Bragg intensities, and calculate an R factor to mea- 
sure its deviation from the measured intensities. Adjust the 
parameters in a direction that will reduce R and repeat the 
simulation; iterate until a good fit is achieved. Following this 
path, the phase solution is unified with the structure fitting; 
it is also natural to unify the i?-factor from diffraction and 
the total energy in a combined objective function to be min- 
imized by the fit. 

Most of this paper is devoted to elaborating the two 
approaches. We begin with the "standard" path in sec. ^ 
focusing on the new "minimum-charge" method, a direct 
method for determining the phases of a general structure. 
This method is not peculiar to quasicrystals, nor to random 
structures. However, the minimum-charge viewpoint is use- 
ful for random structures since it provides an inequality for 
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the diffuse scattering (Sec. [2.2| ) which allows us to quantify 
the degree of disorder without solving the structure, even 
partially. Also, preliminary calculations using the minimum- 
charge approach have suggested the structure of i(AlPdMn) 
is more disordered than any of us had believed. 

As we have just argued, even a perfect solution of the 
phase problem is very far from giving a useful model of the 
structure. This gap can be bridged by the brutally simple 
"factorization approximation": assuming that the intensities 
are those of ideal structure apart from a (Gaussian) "perp 
Debye -Waller" factor. This is critiqued and tested by simu- 
lations in Sec. ^ 

After Sec. |^ we turn to the "unified" path. We first 
review key notions of tile-decorations (Sec. Q), since these 
are used to specify random-tiling structures, practically by 
the definition of a random tiling. Then Sec. ^ describes sev- 
eral fitting procedures, beginning with a recipe for discov- 
ering the appropriate random tiling model by simulations 
(provided a set of interatomic pair potentials is known). Fi- 
nally, we explain in Sec. ^ the special problems posed by 
decagonal structures. 

Our subject is not the long-wavelength "phason" elas- 
ticity which produces diffuse wings around the Bragg peaks. 
These phenomena indeed characterize "random tiling" type 
behavior, but as argued previously in Ref. ljl|], one should 
not necessarily think this accounts for most of the disorder. 
A priori, one expects that the disorder is mostly correlated 
at the scale of a tile edge or so; these correlations are con- 
trolled by the constraints of the tiling rules, and produce dif- 
fuse scattering which is not associated with particular Bragg 
peaks. 



2 Solving the phase problem 

The "standard" approach to structure solution begins by first 
solving the phase problem. This step is somewhat different 
than for ordinary crystals, but it seems not to matter much 
whether the quasicrystal is modeled as ideal or as a random 
tiling. 

We first mention a couple of well-known techniques. 
An old technique takes advantage of the smooth depen- 
dence of the structure factor Fg when plotted as F{g-^), 
where g-'- is the perp-space partner of g. Where |F(g^) 
passes through a zero, one infers that F{g-^) must change 
sign. Another approach, described by Qiu and Jaric, was to 
discover phases using known rational approximant crystals 
(by expressing them as rational cuts of a higher-dimensional 
structure). 

In the rest of this section, we first outline a new ap- 
proach, the "minimum-charge" method (Sec. 2A). Quite in- 
dependent of this, we derive an inequality which is useful as 
a diagnostic of disorder (Sec. 2.2). Finally, we appl y both of 
these ideas to a data set for i(AlPdMn), in Sec. 2.3 , 



imize the average electron charge density. Any periodic or 
quasiperiodic density has the form 



p{r) =Fo+^|Fg|cos(g-r 
where Fq is the average density; |Fg 



I /g and ( 



(2) 



are the 



magnitudes and phases of structure factors. The former are 
obtainable from the measured intensities /g and the latter 



are to be determined. The "minimum charge" methodHll 
considers the reduced density 



~p(v) = p{r) - Fo, 



(3) 



and seeks phases which maximize the minimum value of 
Anin = minp(r). (4) 

r 

The minimum value of Fo which makes p(r) everywhere 
nonnegative is then Fo ~ — Pmin- It is straightforward to 
argue that at an optimal set of N phases (where Fo has 
achieved a local minimum) the minimum pmin is attained 
at 4- 1 exactly degenerate local minima of p{r). Con- 
sequently, for optimal phases, the minimum "background" 
value of the density is found in large regions of the unit 
cell. If some phase assignment allows for the volume oc- 
cupied by this background density to be especially large, as 
we expect when dealing with a true atomic density, the cor- 
responding value of Fq will be especially low. This property 
distinguishes the global minimum of Fo and provides for an 
unambiguous phase solution. 

By sampling p{r) on a "Fourier grid"[[l2|] of a com- 
parable size, the charge minimization problem, say for a 
centrosymmetric space group, is cast into a mixed integer 
programming optimization |]13|] for the unknown (integer) 
signs and (real valued) Fo. Because of linearity (of objective 
function and constraints) there are efficient search strategies 
for this general clas s of problem. For the j(AlPdMn) x-ray 
phases (see Subsec. 2.3 ) the "branch and bound" technique 
was used, where the relaxation of the constraint |sg| = 1 
(on each of the signs) to the weaker statement |sg| < 1 (on 
a subset of the signs) provides bounds on the objective func- 
tion (Fq) that frequently are strong enough to ehminate large 
branches of the search tree. 

There is a similarity between the minimum charge 
(minQ) and the maximum entropy (maxS) methods. Both 
succeed in eliminating or minimizing wiggles in the electron 
density in regions where there should be no charge. The sim- 
ilarity is only superficial, however. On the one hand, maxS 
is fundamentally a strategy for refining phases which are al- 
ready approximately determined, but corrupted by noise; on 
the other hand, minQ makes perfect sense in the context of 
perfect data. Also, minQ was developed specifically for ab 
initio phase determination, whereas the application of maxS 
in this mode has no theoretical basis. 



2.1 "Minimum-charge" approach 

A promising new "direct method" (not specific to quasicrys- 
tals) to solve the X-ray diffraction phase problem was sug- 
gested by one of us[pl|. The key notion of this is to min- 



2.2 InequaUties for X-ray structure factors 

The aim of this subsection is to place a lower bound on 
the amount of disorder, taking advantage of the fact that the 
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density of scattering electrons is nonnegative. This is com- 
pletely independent of the preceding section and, in particu- 
lar, all the conclusions of Subsec. 2.2 can be drawn without 



ever determining or even considering a single phase factor. 
The positivity of p(r) implies the inequality 



|Fg| = |(p(r)e-«--)J<(|p(r)|), = Fo, 



(5) 



where {. . ■)r means a spatial average. It should be empha- 
sized that p{r) already contains an ensemble average: that 
is ju st w hat the Bragg scattering represents, as we noted in 
Sec. |1.1[ Eq. (|]) is the most trivial of the Harker-Kasper in- 
equalities [ p4| . Note, though, that in crystallography, these 
(and later inequalities) were always applied to the solution 
of the phase problem assuming an ideal structure. We apply 
instead as a rigorous diagnostic of a disordered structure, 
using only Bragg (not diffuse) scattering. 

It follows very obviously from (Q) that 
l-fgl/VE < l-Fol/VE l^gP which is conveniently 
written 



< a 



where 



Fo 



(6) 



(7) 



^g#o I-' gl 

Note that a just depends on the first two moments of p, since 



{P)r = 



Fo, 
F2 



El 



(8) 
(9) 



The value of (|g) is that both sides can be related 
to experimentally measurable data, as is evident for the 
Bragg intensities in the left hand side. (The sum omits Fo, 
which is not measurable as a Bragg intensity.) The right 
hand side turns out to depend only on atomic parameters 
such as charges, number density, concentrations, etc., which 
are known from the composition, provided we assume the 
diffracting material consists of a known density of definite 
atomic species (well separated so overlaps of atomic charge 
distributions are negligible). Specifically, suppose there is 
no disorder other than that due to thermal vibrations of the 
atoms. To a good approximation one may then represent the 
electron density as a sum of non-overlapping distributions 
centered around each atom, p^*°"\(r) where i is the atomic 
species. Without knowing the atomic positions, it is possible 
to evaluate the averages in (^ and (^: 



{P)r = n^X^Q^ 

i 



(10) 

(11) 



where n is the number density of atoms, Xi the atomic com- 
position, and 



Q^ 



(12) 



(13) 



and I 



Inserting ( ^0 ) 
%, we conclude 



and (jn]) into the earlier equations, (| 
a = ao, where 



1 



(14) 



- 1 



is expressed purely in terms of the atomic parameters n, Xi, 
Qi and P^, 

Our inequality (^ is rigorous. If it is violated by the 
value of ao computed from the data, it means we were mis- 
taken in our assumption that the electron density p(r) has 
100% occupation of every site and no chemical disorder. Let 
us repeat the derivation above, assuming a distribution of oc- 
cupations (in the case of a single species for simplicity). In 
place of dl4|) we find 



(15) 



in the case where a <C 1. Here z{y) is the distribution of 
the fraction of charge from sites with fractional occupation 
y. Disorder decreases {y) and thereby, according to (15), in- 
creases a beyond the value in (|l4). If we allow for enough 
occupational disorder, inequahty ( Sb will be satisfied. 



2.3 Example: the i(AlPdMn) density map 

Recently [|l5|] the "minimum-charge" method was applied to 
determine the phases for i(AlPdMn), which has the cen- 
trosymmetric space group F53^, from the x-ray data of 
Boudard et al.JT^. This data set has 503 symmetry inequiv- 
alent reflections, 360 having (Tg//g < 0.5; we used all the 
reflections. 

Figures |l] and ^ both show the reconstructed elec- 
tron density using signs obtained with the minimum charge 
method. Slices through physical space (Fig. |l]), and along 
rational planes (Fig. both show disorder, as we argue in 
the next two paragraphs. Disorder of this sort is expected in 
a random tiling, but our observations here do not rule out 
other origins. 

Figure |l](a) shows the electron density in a 2-fold 
(mirror) plane of physical space, centered on the "no" 
atomic surface. While a subset of the atoms corresponds ex- 
actly to the Mackay cluster second shell (12-1-30 atoms), we 
see that there are also several pairs of well defined atoms 
having unphysically short separations. Similar evidence of 
split positions has been noted before by others, but was 
never ruled definitive because of the possibility of being a 
truncation artifact|]l7|] . Because of the flatness of the back- 
ground in these reconstructed densities we feel that trunca- 
tion cannot be invoked to explain these examples of split 
positions. A preliminary survey of the reconstructed density 
suggests that split positions are actually quite common. 

In the periodic 5-fold plane (Fig. ^ we see three 
atomic surfaces centered at Wyckoff positions with icosa- 
hedral symmetry. The positions and net charge in these sur- 
faces are consistent with the "spherical model" of Boudard 
et q/JT^. We observe, however, that the density profiles 
of all three surfaces are very rounded, rather than step- 
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function-like. |^ The width of the appropriate broadening 
function in the factorization approximation (see eq. ([T^, 
below) would be comparable in magnitude to the diameter 
of the largest (ideal) atomic surface in the spherical model. 
This explains, for example, why the sizes of the largest 
("no") and smallest ("6ci") surfaces of Ref. [ p^ are com- 
parable in our image, whereas in the spherical model their 
diameters are in the ratio 2.3:1. 

Even before determining the phases, the Boudard 
datajl^] can be manipulated to show (using the inequali- 
ties in Sec. ^!2[ ) that the i(AlPdMn) phase possesses con- 
siderable partial/mixed occupational disorder. To evaluate a 
from eq. (p^, we first approximated the atomic form factor 
for each species z as a Gaussian, 



(16) 



with the parameters Bi chosen to reproduce the integrals 
Pi in ( p3[ ) obtained using Hartree-Fock wavefunctions. This 
gave Bai = 0.012, Bpd = 0.0042, and Bm„ = 0.007 (in 
units of R^). Since Bpd is close to the measured Debye- 
Waller factor for Pd[|l8|], Bdw = 0.0044, the thermally av- 
eraged distribution for Pd is better represented by a Gaussian 
with B = 0.0086. Since Pd atoms make the main contribu- 
tion to a, we used the same Debye- Waller factor to correct 
the Al and Mn distributions. Using the measured composi- 
tion and density [pq], we finally arrive at ao — 0.0439. 

The sum ) ] [fgP in (^ is in principle determined by 
a measurement of all the Bragg intensities. Barring surprises 
in the unexplored regions of reciprocal space, we believe the 
intensities that have been measured exhaust this sum: the 
weaker measured peaks of our data set already make a negli- 
gible contribution to it, and peaks with larger g or are cut 
off by Debye- Waller factors. Our conclusion about the dis- 
order follows from the fact If we assumed no disorder, then 
OL — and the three most intense reflections (18/29,52/84, 
and 20/32) would violate inequality (||) by as much as 40%. 
Our estimate of disorder is certainly conservative since there 
is no reason to e xpect that inequality would be close to 
an equality. (Sec. 3.1.1 also indicates a large disorder in this 
material.) 



3 Factorization approximation 

The determination of the phases, and hence of the averaged 
scattering density, is not the end of the structure determi- 
nation in the random-tiling case. A process of simulation 
and fitting is still called for. (See Sec. |[ below.) The fac- 
torization approximation - equivalently the "Phason Debye- 
Waller factor" - is a shortcut of dubious validity, which is 
nevertheless attractive since it offers the hope of relating the 
data to an ideal model with much less effort. 



We caution that some important details are lost when the weak 
reflections are omitted; in other words, truncation error may introduce a 
spurious apparent randomness. For example, with onlv 300 reflections 
the profile of the atomic surface at the center of Fig. WO) looks quite 
Gaussian. (Its correct appearance is two-peaked because its center in 
perp space is occupied by Mn, while the rest of the atomic surface is Pd 
which possesse twice as many electrons as Mn.) 



In fitting diffraction, it was natural to generalize the 
Debye-Waller factor to perp space and to assume 



^rand _ ^idoalp — W(g) 



(17) 



for the Fourier amplitudes of the random and ideal tilings 
(labeled by "rand" and "ideal" henceforth). Here the "perp 
Debye-Waller factor" is given by 



iy(g) 



72 



(18) 



where g^ is the perp-space wavevector corresponding to g. 
We will call this the "factorization approximation". In direct 
(perp) space, in view of (|]), eq. (|7|) is equivalent, for an 
icosahedral quasicrystal, to 



P 



rand 



(19) 



where w(u^) is a (normalized) smearing function with a 
Gaussian profile: 



Zi;(u''^) = (27ro- 



2N-3/2 



cxp[-|u^|V2a2] 



(20) 



Notice that it follows from ( [l9| ) that the random tiling ver- 
tices have a perp variance increased over the ideal tiling by 



(|h^P},and-(|h^|')idcal=3a2 



(21) 



For a general quasicrystal, "3" must be replaced by "dj^" 
in eqs. ( [l9| ) - (|I]). (We wrote them assuming perp space is 
three-dimensional as in the icosahedral case.) 

It must be borne in mind, however, that there is no 
exact basis for the factorization approximation. One may be 
misled because elastic theory tells us all long wavelength 
perp fluctuations are Gaussian. But this relation is not use- 
ful, since we expect much of the random-tiling disorder to 
be short-range. 

Another way one may be misled is that, when the 
higher-dimensional cut construction is overemphasized, a 
natural cartoon of the random tiling is to take the same 
higher-dimensional crystal as in the ideal case, but to per- 
mit the cut-surface to undulate - to deviate from being a flat 
plane. If the cut-surface has Gaussian fluctuations that are 
completely uncorrelated with the atomic surfaces, then the 
factorization approximation would become exact. This "un- 
dulating cut" notion was critiqued in Sec. 4 of Ref. [jl[].0 

The factorization assumption is fundamentally ill- 
defined since the same random tiling (e.g., of Penrose 
rhombi) may be obtained as the randomization of various 
different quasiperiodic tilings (e.g. the generalized Penrose 
tilings in two dimensions) - yet eq. ( p^ ) demands that we 
identify a particular ideal tiling as the one which has been 



^ This does imply a formula like (|l^ relating p''''"'*(h-^) to the 
coarse-grained distribution p^*^""^* (h^), obtained after averaging the 
"center-of-mass" in perp space locally, from a region several tiles wide, 
rather than over the entire system. 

^ In places where a straight cut would always cross exactly one out 
of two atomic surfaces, an undulating cut may cross both or neither of 
them, which produces an overlap or gap between parts of tiles. Addi- 
tionally, two slightly different undulating cuts often produce the same 
real structure, so one must be careful in assigning a statistical weighting 
to the ensemble of undulating cut surfaces. 
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randomized. The physical meaning of ( [T^ ) is dubious, too: it 
says that, when the tiling is randomized, the vertices with a 
prescribed local environment, which in the ideal tiling have 
a specific value, now get shifted by a random offset 
which has a probability distribution (pO|). But there is no rea- 
son the random variable should have the same variance 
( [21] ), independent of what kind of local environment is in 
question or of how large is its location. 

We claim the factorization approximation works triv- 
ially for sufficiently small g^. After all, an integral of form 
^ can always be written as exp(— Co + |C2|g^P + ••■)' 
where C2k is the cumulant of order 2fc of the density p(h^). 
Then since the 2fc = 4 term in the exponent is of order 
(g-'-)'', it follows that F(g-'-) has an roughly Gaussian shape 
at small g^, even for the ideal tiling (see Ref. [^). The real 
test of the factorization approximation is whether it works 
for larger g^, with the fitted in ( 18) being independent of 
g^. Surprisingly, this seems to be true for our simulations, 
as the rest of this section shows. 

3.1 Simulation tests 

We carried out simulations of the completely random tiling 
of rhombohedra [ pO[ pl| , p2[ ], as a toy model to test (for the 
first time) the relation between diffraction amplitudes in the 
ideal and random tilings: in particular, to test the validity 
of Eq. (p7|). These measurements are preliminary; we hope 
they can be repeated in the future, more systematically and 
perhaps on models closer to the real quasicrystals. The rest 
of this section is devoted to these simulations. 

We used the cell of the "8/5" cubic approximant with 
periodic boundary conditions, which contains 10 336 ver- 
tices (an equal number of rhombohedra). For equilibration 
we allowed 500 Monte Carlo steps (MCS) / vertex (these 
are flip attempts of which only 17% cause flips), and then 
our ensemble for diffraction consisted of 1000 sample taken 
at intervals of 250 MCS/vertex. (This is adequate to decor- 
relate this maximally random model, ^ 

but quite inadequate when there is a Hamiltonian 

iHi.) 

A novel detail of our simulation is that periodic 
boundary conditions are also assumed, rather arbitrarily, in 
the perp direction. Thus, only a finite set of (178)"^ val- 
ues is possible. In practice all vertices fall within a smaller 
domain of these in perp-space - ISO'^ points were used here 
- and configurations are efficiently represented as a lattice 
gas on this grid. 

3.1.1 Results on perp-space distribution 

It is possible to visualize the smearing of eq. (|l9|) directly, 
rather than through the perp Debye- Waller factor. On the 
theory side, it is straightforward to construct a histogram 
of values for the vertices of the random rhombohe- 



dral tiling (see Sec. 3.1, below). A random-tiling deco- 
ration model [Bsl] has been proposed for i(AlCuFe) and 



^ The relaxation time of the slowest Fourier mode in the 8/5 cell was 
found to be only 280 MCS/site [Q. Our Fig. ^ agrees well with Fig. 1 
of Ref. |§; furthermore, a different run with the same protocol yielded 
elastic constants {K-i , K2) = (0.84, 0.50) within 2% of the results for 
the same size in Ref. Fig. 2. 



i(AlPdMn). This is based on the rhombohedral random 
tiling, which has simple-icosahedral space group. However, 
the atomic decoration treats even and odd vertices inequiv- 
alently, so the atomic structure is face-centered icosahedral, 
like the real quasicrystals, and is similar to prior models 
based on diffraction [^, |l^ In particular, the bci atomic sur- 
face is (in this model) made up precisely of atoms that deco- 
rate the even rhombohedron nodes and should have the same 
perp-space distribution /o(h^) as all nodes do in the simu- 
lation (since even and odd nodes have the same pih.^) in 
this model). In short, the experimental perp-space distribu- 
tion of the hci surface is just a section through the density 
map (Sec. ^. 

Surprisingly, assuming the phase reconstruction is 
correct, i(AlPdMn) - the most perfectly ordered quasicrys- 
tal to date - appears "more random than random" as quanti- 
fied by the apparent value of in Fig. ^ This conclusion is 
also supported by the quite small numerical values of elastic 
constants in i(AlPdMn), as measured from diffuse scatter- 
ing JT^, in comparison to models [p2|]. We know only two 
possible explanations of this disorder, neither of which is 
very convincing: 

Explanation ( i): The true tiles are deflated compared 
to the assumed ones by some factor like (such 
rescalings are discussed in Sec. 9.2.1 of Ref. [|I]]). But 
for i(AlPdMn) this seems impossible, as the true tile 
edge would be no larger than an interatomic spacing. 
Explanation ( ii): The ratio of elastic constants is close 
to K2/K1 — —0.75 at which a "phason instability" 
occurs, with a divergence of perp fluctuations Pq]. 
Within the entropic-stabilization scenario, such an in- 
stability is indeed expected when the temperature is 
lowered [j^ - but only near a critical temperature, not 
over a wide range. 

Further understanding awaits a better understanding 
of the currently unsettled experimental situation. In partic- 
ular, Ref. found that the diffuse scattering differed by 
a factor of 20 between two slightly different samples of 
i(AlPdMn); it is unclear which of these is hke the sample 
of Ref. @. 

3.1.2 Computing Bragg and diffuse scattering 

The rest of our studies in this section depend on the diffrac- 
tion. To compute this, each vertex was taken to be a point 
scatterer with form factor unity. In the size 8/5 approximant, 
computing the Fourier transform at every wavevector in the 
interesting range would take weeks or years of CPU time, so 
we used lists of selected wavevectors (see below). 

When one uses periodic boundary conditions, of 
course the real reciprocal-space vectors fall on a discrete 
grid. Due to our special adoption (see above) of perp peri- 
odic boundary conditions, the perp reciprocal-space vectors 
also lie on a discrete grid. Indeed each of our real reciprocal- 
space vectors q is the projection of infinitely many 6D 
reciprocal-lattice vectors, and can thus be considered as a 
Bragg vector (with the appropriate distortion since this is an 
approximant). We select the smallest of these to make the 
correspondence unique. In practice most of the q vectors 
correspond to such a large g^ that the Bragg amplitude is 
zero to computer precision. 
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We defined the "Bragg" and "diffuse" parts of the 
diffraction according to their behavior in the infinite-size 
limit, so we need a new operational definition or procedure 
to separate the Bragg from the diffuse intensity in a fixed, 
finite lattice. This is easy to do, in view of the discussion 
in Sec. LI: the amplitude's ensemble average is taken to be 
the Bragg amplitude, while its variance is the diffuse inten- 
sity. However, even if each sample were very similar to the 
previous one, it is offset by a random fraction of the sim- 
ulation cell, independent of the previous offset. Hence, be- 
fore taking the Fourier sum over each sample, we located 
the "center-of-mass" of its perp-space coordinates, and then 
uniformly shifted the vertex coodinates (by integer multiples 
of the underlying grid of the lattice-gas) so as to make this 
center-of-mass nearly zero. 

In this fashion we could measure Bragg intensity even 
on peaks where the diffuse intensity was perhaps ^ 10 times 
as big. ^ The results (along a twofold axis) are shown in 
Fig. Bragg peaks of the ideal tiling are identifiable down 
to amplitude 10~^ relative to the maximum (far smaller 
than experimentally measurable). Bragg amplitudes of order 
10~^ in the ideal tiling became of order 10^"' in the random 
tiling, which is about where they become lost in the diffuse 
noise. In the best experimental data sets, the smallest mea- 
sured Bragg amplitude is down about 10^^ from the largest 
one. 



3.2 Results: perp-space Debye- Waller factor 

Fig. ^ shows the Bragg intensities plotted against |g^| for 
selected wavevectors; the amplitudes (same data, but with 
signs) are plotted in the upper panels of Fig. |[ We can define 
the ratio 

which is plotted against Ig"*"]^ in Fig. ^(a). The straight-line 
fit confirms (|l8|) with ct^ = 0.145. On the other hand, the 
variance (|h^f ) has been (easily) measured for both ran- 
dom and ideal tilings: inserting these numerical values into 
Eq. (|l|) gives 3^2(1.670) - (1.236) = 3(0.148), in per- 
fect agreement. The only caution to be mentioned is that the 
deviations which appear small because of the logarithmic 
scale of Fig. ^ may in fact be far greater than the measure- 
ment errors of the intensities. By the time the factorization 
approximation breaks down completely, at jg"*"]^ ~ 25, the 
random-tiling Bragg peaks (see top panel) are too small to 
measure in a real experiment. 

Do the Bragg amplitudes behave differently when 
there are several orbits with different species of atoms, or 
when the atoms do not sit on vertices? To check this, we also 
calculated the diffraction from a decorated model in which 
atoms with form factor + 1 sit on all rhombohedron vertices 
and another species with —1 sits on mid-edges and also on 
prolate rhombohedron axes at ideal points corresponding to 



body centers of six-dimensional cubes. |] This places 54 120 
atoms in the "8/5" size simulation cell. 

The answer - appearing in Figures ||, ^ and ^- is that 
the decorated model is the same qualitatively - for example, 
is fit to 0.143, essentially the same value. (The fit is just 
a little bit worse.) 

3.2.1 Correctness of the phase? 

A special concern was whether the random-tiling amplitude 
always has the same sign as the ideal-tiling. Recall that, as a 
function of |g-'- 1, the ideal-tiling amplitude oscillates and has 
zero crossings at certain wavevectors. The first zero-crossing 
of the random-tiling vertices agrees precisely with that of the 
ideal-tiling vertices; however, subsequent crossings at larger 
|g^| are noticeably different. This indicates that the perp- 
space profile of the atomic surface is not exactly given by 



3.3 Integration over diffuse wings 

Some of the intensity lost to diffuse scattering is recov- 
ered when one integrates over the diffuse wings around each 
Bragg peak; some such integration is always performed due 
to the resolution function. We did another simulation to test 
whether this, in some way, undoes the DW reduction. 

We summed up the total intensity in a sphere centered 
on each Bragg point. The sphere radius was 3 grid spac- 
ings of the lattice of allowed wavevectors, so it contains 123 
wavevectors. (The spacing of our grid of discrete q points is 
reciprocal to the cubic lattice constant 18.87 of the 8/5 ap- 
proximant. i.e. |Aq| = 0.333 in units of inverse rhombohe- 
dron edge. As is visible in Fig. Q the separation between no- 
ticeable Bragg peaks is 6 or 10 grid spacings, i.e. t^^tt or tt, 
so these spheres do not overlap.) Due to computer time lim- 
itations, we only performed this computation for 27 Bragg 
vectors, selected to be representative of the range of val- 
ues - a very small fraction of all the usable Bragg vectors. 

Fig. ^ also shows these integrated results (crosses). 
Unexpectedly, the integration seems to restore the same frac- 
tion c of the lost intensity of every peak independent of g-'- 
- the same when |g^ | is so small that the loss is invisible on 
the figure, and out to peaks where the ideal intensity is down 
by ^ 10"^ and the random Bragg intensity is down a further 
^ 10^^ from the ideal. (Obviously, though, c depends on the 
radius of our integration sphere.) The formula for this is 



c/sadoal + (1 - c)Ib,-. 



(23) 



where /s. ideal and /strand are the Bragg intensities for the 
ideal and random tilings, respectively, and //,iand is the in- 
tegrated intensity for the random tiling. To get a good fit, we 
had to include the constant D which represents some diffuse 
intensity spread uniformly over reciprocal space. |^ Formula 



If the system size were changed, then of course the Bragg amplitude 
and the diffuse intensity both scale linearly with the system's volume. 



^ The positions are a simplification of Ref. |p4|, but the assignment 
of species is different. Negative form factors are realized by neutron 
diffraction for some isotopes. We chose these form factors to make this 
diffraction pattern as different as possible from that of just the tile ver- 
tices. 

* Uniform diffuse intensity implies, of course, a component of the 
density that is completely uncorrelated from point to point. We know 
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( p3| ) is plotted for comparison in Fig. q|; Fig. ^confirms more 
directly its remarkable success (which is somewhat worse in 
the decorated model). 

It is intriguing to speculate how a fit using the inte- 
grated intensities would behave if (^) were exactly correct 
for real experimental data, particularly if (due, say, to ex- 
tinction problems) the strong peaks at small were omit- 
ted from the fit. Then /s,tand could be neglected in ( p3[ ) and 
^/,iand /s, ideal, i-C-, the ideal structure would be recon- 
structed (apart from a wrong factor of c in the absolute nor- 
malization of the density). 

Note that, inserting the phason Debye -Waller hypoth- 
esis of (|T^ and (|8|) into (g) impUes 

//a-and//B,idcal « C + (1 - c) (exp -a^\g^\^) (24) 

Thus the integrated intensities are not expected to be fit- 
ted by the phason Debye-Waller form //.rand/^s, ideal ~ 
exp(— cr||g-'^P). However, if the |g^| values are so small 
that exp(— cr^lg-'-p) « 1 — cr^ |g^p, the fit will give an 
apparent cr| « {1 — c)a^. Indeed, a smaller value was 
fitted by de Boissieu for low-resolution data than for high- 
resolution data in Ref. [^. 

3.3.1 Overall transfer to diffuse intensity 

A key issue in Sec. ^ was what fraction of intensity is trans- 
ferred to diffuse scattering in a random tiling. We investi- 
gated this using the same random rhombohedral tiling. We 
summed the total diffuse intensity and the total intensity for 
all wavevectors lying in the 2-fold plane. The fraction of 
the total intensity that was "diffuse" was 0.145 in the 3/2 
approximant and 0.095 in the 8/5 approximant: the 3/2 ap- 
proximant shows a strong finite-size effect. This fraction is 
artificially low in the 2-fold plane which contains many re- 
ciprocal lattice vectors with small g^. The fraction of all 
intensity which was diffuse was 0.388 in the 3/2 approxi- 
mant; it was too much to compute in the 8/5 approximant, 
so we can merely guess it is ^ 0.25 (based on the trend in 
the 2-fold plane). 

3.3.2 Digression on diffuse wings 

A great deal more can be learned if one measures the de- 
tailed shape of the diffuse wings at high resolution, in- 
stead of just integrating over them. Diffuse scattering, un- 
like Bragg scattering, depends on the spatial correlations of 
disorder. 

Elastic theory governs the long-wavelength fluctua- 
tions of (r) in random-tiling-like quasicrystals; energy- 
stabilized models seem not to predict a gradient-squared 
elasticity . Also, the shape of the diffuse wings of Bragg 
peaks is a key prediction of the elastic theory l|2^, . 
Thus, the match between experiment and elastic theory in 
i-AlPdMn and in i-AlCuFe, [ |3l| ] is evidence that these qua- 
sicrystals have a random-tiling nature. (Some of us recently 
commented on these experiments and calculated some elas- 
tic constants, in p^ ] and [p2[].) 



no specific defect tiiat does this in our random tiling; the alternation 
between the two ways of packing four rhombohedra into a rhombic 
dodecahedron might approximate this. 



However, such analysis of diffuse scattering is out- 
side the scope of this paper, which is crystallographic inves- 
tigations based purely on Bragg intensities of many peaks, 
aimed at discovering the microscopic atomic arrangements. 

3.3.3 Global diffuse scattering? 

An unresolved question, and of relevance to background 
subtraction, is the diffuse background far away from Bragg 
peaks. It has been suggested by Ref. | p7[ ] that one can ac- 
count for all of this by summing over the diffuse wings of 
all Bragg peaks, as approximated by elastic theory | p6[ p9|]), 
even from those that are farther away in reciprocal space 
than the typical separation of strong Bragg peaks. This is 
surprising since the elastic theory is a long-wavelength, 
coarse-grained theory that ought to break down at the cor- 
responding real-space distances (of the order of a tile edge). 

In any case, this formula cannot be integrated over all 
of reciprocal space, since even a single diffuse wing has the 
schematic form J d'^q(l/iir|qp). It is of interest to write 

(|h I ) = i2n) d q (25) 

where the 3 comes from taking the trace of a 3 x 3 ma- 
trix, and e{K2/Ki) is a dimensionless factor which aver- 
ages out the angular dependences. For K^jKi ^ 1, it can 
be shown that g 1 + {%l^){K2l K^f + . . ., while it di- 
verges for K2/K1 — » 3/4 (the usual modulation instabil- 
ity [Hi). Noting iha\_{Ki,K2) = (0.8,0.5) for the rhom- 
bohedron tiling | pO[ , pl| and e « 2, we obtain Qmax ~ 3.5 
which is comparable to tt (the spacing between rather strong 
Bragg peaks). 



4 Review of tile decoration models 

Before going on to the "unified" approach, we will need 
some more terms and concepts - associated with tilings and 
their decorations (in real space) - to be used in the two sec- 
tions following. To start off, we imagine that an ensemble 
of certain "low-energy" atomic configurations adequately 
represents the states of the quasicrystal at temperatures be- 
low lOOOK. Then we insist that every "low-energy" con- 
figuration must correspond to a unique tiling. (Hence we 
sometimes say "tiling" as a shorthand for "atomic config- 
uration represented by the tiling," e.g., when we speak of 
"the tiling's energy.") 

4.1 Decoration 

The mapping from the tiling to an atomic configuration is 
called a decoration, since each tile of a given class has atoms 
placed in the same sites on the tile. If necessary, this place- 
ment is allowed to have context dependence on the neigh- 
boring tiles. But if the site is sometimes free to receive either 
of two species, depending on neighboring tiles, it may be 
preferable to consider this tile as having two "flavors" each 
with a unique decoration. Sometimes atoms are assigned as 
decorations of other "tile-objects" such as vertices, edges, 
or faces, in order to reduce the number of site classes and 
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numerical parameters (see Sec. II A of Ref. |]33[]). In all this, 
the principle is that all the constraints and correlations are 
ascribed to the tiling level, so that the statistics of the tilings 
contain implicitly all the statistics of the atom configura- 
tions. ^ 

In random tilings, a rearrangement which turns one 
valid tile configuration to another, while affecting only a 
small cluster of tiles, is called a tile flip. For a good deco- 
ration rule, such that all tilings in the ensemble have sim- 
ilar energies, there is a substantial agreement between the 
atomic configurations on those tiles before and after the tile 
flip. (If this weren't so, one of the two configurations would 
have a much higher energy due to interactions with atoms in 
neighboring unflipped tiles.) Often, only a couple of atoms 
change positions in a tile flip. 



4.2 Supertiles 

The close overlap between the "before-flip" and "after-flip" 
atom configurations often gives us freedom to change the 
size of tile while only slightly changing the structure model. 
When switching from large tiles to smaller tiles, the price 
of using small tiles is some simplifications in the decoration 
rule, and/or additional constraints among the tiles; the price 
of using large tiles is having many sites on the tiles with in- 
dependent parameters. Replacing tiles by supertiles always 
entails a rebinding of their decorations. 

The supertiling phenomenon is particularly prominent 
in decagonals and is already familiar experimentally - par- 
ticularly in the various subphases of (i(AlNiCo) [ p4| ] - and 
is also known as a cause of mis-indexingl|3^]. The "inflated" 
tiles in Penrose's tiling are a special case of supertile, but 
most supertiles are not literal inflations. Physically, starting 
from small tiles with a tile Hamiltonian, the small energy 
differences between tilings favor a sub-ensemble of degen- 
erate lowest-energy configurations and often these may be 
represented as tilings with bigger "supertiles", each deco- 
rated in a single way with the original tiles; a special case is 
the interactions which favor a local pattern or "cluster" 
Further remarks on supertilings are in Ref. ||37 



4.3 Recipe to discover basic structure 

To proceed, one must first decide which tiling geometry best 
repr esents the structure and its degrees of freedom. Table 
4.3 lists a menu of available random tilings which have been 
used in decoration models of real materials. (Some of the 
decorations were designed for a quasiperiodic tiling, but are 
compatible with the random one.) For the tilings marked 
"no" in column 3 of the table any simulation is just barely 
practical, since each update move must involve an indefinite 
number of tiles, [fl H ISOll . 



Decoration models have been formulated such that decorations of 
different tiles sometimes produce two atoms in the same position, the 
rule being that these get replaced by one atom. We do not consider such 
decorations, because in analyzing or optimizing sums over the struc- 
ture, giving either diffraction amplitudes or total energies, it is awkward 
when one cannot assign each atom to a well-defined tile-object that it 
decorates. 



Table 1. Important tilings for realistic structure models. Under "Lo- 
cal?", a "yes" means a local update exists; only such tilings are easy to 
simulate. (We have sometimes called the hexagon-boat-star tiling "Two- 
level", or the rhombohedral tiling "3DPT"; the latter may also contain 
rhombic dodecahedra.) Under "symmetry", ico, 10, and 12 are short for 
icosahedral, decagonal, and dodecagonal, respectively. References are 
indicated for tilings and decorations. 



Tiling 



Symm. Local? decorations 



rhombohedral [ [20[ [Zl|1 
canonical cells p^^o[ ^l| 
binary MSl M] 
hexagoi>boat-star [ttffl 



ICO 

ico 
10 
10 



yes 
no 
yes 
yes 



):| 



rectangle-triangle [|50p 10 
square-triangle tesl 12 



j(AlPdMn. 
i(AlMn): [p 
d(AlCuCo 
d(AlMn): 
d(AlPdMn 
d(AlCuCo): 
d(AlNiCo) 
d(AlPdMn): 
d(AlNiCo), 
dd(VNiSi) ' 
dd(Tai.6Te 



The proper size of the tiles is not self-evident even if, 
say, one already knows the exact structure of a large approx- 
imant crystal. Consider the most extreme case, a structure 
that looks just like the Penrose tiling. That structure may 
be broken into tiles at any level of inflation; furthermore, 
at each level, the tiles may be represented as Fat/Skinny 
rhombi, as Kites/Darts, as Hexagons/Boats/Stars, or in a 
couple of other ways. At which level can we break apart the 
tiles and put them together differently, to build other approx- 
imants or random tilings? Only a computation (or intuition) 
of the energies can tell us. 

In the rest of this section, we outline a recipe to find 
the right tiling, and a coarse version of the structure model, 
before using any diffraction data, from simulations of lattice 
gases [p^, or of random tilings with small tiles, when mi- 
croscopic Hamiltonians are available. There are four or five 
stages. 

4.3.1 Total energy code 

Before starting, one must have a means to compute a total 
energy for any plausible configuration of the atoms forming 
the quasicrystal. Normally, the best choice is an effective 
pair potential fitted to, or derived from, ab-initio data [^^. 
In principle one can use ab-initio energies directly [Q, but 
these are limited to very small system sizes. 

4.3.2 Lattice gas simulation 

Construct a list of (available) sites, forming a quasilattice. 
The quasilattice constant should be smaller than the inter- 
atomic spacing, so that the occupied fraction of sites is fairly 
small compared to unity. A Monte Carlo simulation is car- 
ried out (using the Hamiltonian described in Stage 1), in 
which atoms are allowed to hop as a lattice-gas among these 
discrete sites; (usually) the temperature is gradually reduced 
to zero and the final configuration is examined by eye. p7|]. 
This annealing must be repeated over and over, since the 
system will end up in different final configurations. (That is 
natural when there are many nearly degenerate states, sepa- 
rated by barriers.) 

Although the states are not constrained to be tile deco- 
rations, it may turn out that the low-energy states can be rep- 
resented as such. Discovering this representation is an art. 
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and there may be more than one correct answer (in that dif- 
ferent tiling/decoration combinations might generate iden- 
tical atomic structures.) The set of possible local patterns 
found in the resulting tilings implicitly defines a packing or 
matching rule for the tilings. Thus, the lattice gas simulation 
serves as a systematic procedure to discover both the "right" 
tiling (with its rules) and the "right" decoration (relating the 
tiling to atomic positions). 

A hybrid random tiling-lattice gas simulation signifi- 
cantly improves on the plain lattice-gas simulation. The sys- 
tem has two kinds of degrees of freedom (both discrete). The 
first kind is a random tiling (which does not need to be the 
same as the random tiling we are trying to discover by this 
technique.) Each tile is decorated in a deterministic way by 
candidate sites: these sites take the place of the quasilattice 
defined in the plain lattice-gas technique. The second kind of 
degree of freedom is a lattice gas on these sites. The Monte 
Carlo simulation must allow tile flips as well as atom hops. 

A different kind of hybrid simulation is to use a ran- 
dom tiling of rather small tiles, with a deterministic decora- 
tion by atoms. In effect, this is similar to a lattice-gas sim- 
ulation, in that (with around one atom per tile), practically 
every topological arrangement of atoms is represented by 
some tiling. 

4.3.3 Simulations of decorated random tilings 

At this stage, the degrees of freedom are tiles, and each 
tiling is decorated deterministically with atoms. A Monte 
Carlo annealing is perfomed in which all tilings are allowed, 
but are weighted as usual by the interaction energy of their 
atoms. Different variations of the decoration rule are tried 
out, in an attempt to find the rule which gives final states 
with the lowest energy. This stage, like stage 2, refines dis- 
crete degrees of freedom. Also like stage 2, it finds both a 
decoration rule and tile packing rules, but in a less tentative 
fashion. 

Having carried out stage 3 in small approximant crys- 
tals, one then repeats it in larger approximant crystals. (Data 
from the small approximants could be misleading as there 
are many local patterns of tiles which can fit in a large ap- 
proximant or an infinite tiling, but not in a small approxi- 
mant.) 

4.3.4 Optimization of continuous parameters 

In this final stage, the decoration is fixed, but decoration pa- 
rameters such as the atomic coordinates on the tiles are ad- 
justed so as to minimize the total energy. One should still 
decorate more than one of the final annealed tilings from 
stage 3. 

After stage 4, the decoration and its parameters are fi- 
nally fixed; the only degree of freedom is tile configurations. 

4.3.5 Construction of tile Hamiltonian 

The energy can be calculated for each possible tile configu- 
ration by decorating it with atoms and using the interatomic 
pair potentials, but this is time-consuming for a large-size 
approximant. It may be possible to find, and use, a "tile- 
Hamiltonian" which is an explicit (and simple) function of 
the tiles or pairs of neighborin g til es, which accurately mim- 
ics the pair-potential energy [ p3| , |4^ . This stage could be 
skipped. 



5 Unified fit from simulations 

The "unified" approach is the second of the two paths for 
structure determination described in this article. It may be 
applied after one has a rough or moderately detailed picture 
of the atomic structure and of the tile degrees of freedom 
(from related approximants, for example). As noted in the 
first path, the procedure for fitting the averaged density, af- 
ter the phase problem is solved, involves fits very similar to 
those required by the "unified" method, and thus does not 
deserve a separate discussion. 

5.1 Random tiling simulation 

An honest approach to a random tiling - more precise than 
the "factorization" approximation of Sec. ^- requires a sim- 
ulation, which can be used in several different modes. Mode 
(i) would be to assume a priori a particular random tiling 
ensemble (that means a fixed list of configurations sampled 
from a simulation); only the decoration parameters are var- 
ied. Mode (ii) would be to hold the decoration parameters 
fixed and vary the "tile Hamiltonian" parameters to con- 
trol tile correlations, until an optimum fit is achieved. In 
this case, every iteration requires a fully equilibrated Monte 
Carlo simulation of the random tihng; the approximant cell 
size may be limited by time constraints since we must do 
repeated Monte Carlo simulations. Mode (iii) would be a 
"simulated annealing" of the tiling: one performs a Monte 
Carlo simulation in which the i?-factor itself replaces the 
tile Hamiltonian in the Boltzmann factor. As the tempera- 
ture is reduced, only the tile configurations which optimize 
the i?-factor will be represented. 

5.1.1 Structure determination assisted by energy 
calculations 

We strongly believe in the need to combine diffraction and 
energy inputs to structure determination, because they con- 
tain complementary information. For example, some rare 
atoms are impossible to pin down by diffraction, but very 
easy to decide on the basis of energies. Again, transition- 
metal atoms from the same row ordinarily cannot be distin- 
guished by X-ray diffration, but they can be clearly assigned 
by the use of potentials [^^. 

In fact, exactly the same strategies can be applied ei- 
ther to the energy values or to the diffraction data. For ex- 
ample, one could run through all the steps of the "Recipe" in 
Subsec. 1-.3 minimizing, instead of the total energy, the R- 
factor for the fit to a diffraction data set. (In a diffraction fit, 
the parameters appearing in the Stage 3 description might in- 
clude Debye -Waller factors or partial occupations.) While it 
is desirable to run energy and diffraction fits separately, the 
"most realistic" fit should do them in parallel, minimizing 
some linear combination of the total energy and the i?-factor 
(or the x^) quantifying the mismatch to the measured Bragg 
intensities. (This would be similar in spirit to the method of 
"least-squares with energy minimization" in ordinary crys- 
tallography 



^0 The familiar invocation of "steric constraints" to rule out unphysi- 
cally close atoms could be considered an informal application of energy 
information. 
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5.2 Technical issues in simulations for diffraction 
fitting 

Here we discuss how one might, in practice, do an iterative 
calculation which combines (i) Monte Carlo reshuffling of 
the random tiling; (ii) computation of Fourier sums to obtain 
the structure factors; and (iii) adjustment of parameters to 
optimize the fit. 

5.2.1 Simulation cell as an approximant 

A simulation, of course, can only be done in a finite sys- 
tem, normally with periodic boundary conditions. Period- 
icity forces a net background phason strain; this is mini- 
mized when the simulation cell has the size and shape of 
the unit cell of a good rational approximant of the quasiperi- 
odic tiling. To calculate the diffraction, we repeat this cell's 
configuration throughout space. Hence reciprocal space has 
a very fine grid of Bragg peaks, but most of them have 
negligible intensity even for an approximant of an ideal 
quasiperiodic tiling. Of the Bragg peaks identifiable in the 
ideal case, most are lost in the random-tiling case. (Their 
intensity is greatly reduced by what is loosely called the 
perp-space Debye- Waller factor, while meanwhile a uniform 
diffuse background appears almost everywhere in recipro- 
cal space; the Bragg intensity is overwhelmed if it is several 
times smaller than the diffuse intensity.) 

Each orbit of equivalent quasicrystal Bragg peaks 
breaks up into several orbits of inequivalent Bragg peaks 
of the approximant. We should average their amplitudes, to 
reduce the systematic error (due to phason strain) and the 
statistical error (from the finite size or time of the random- 
tiling simulation). Before doing such an averaging, one must 
be careful to shift the unit cell so that the different approx- 
imant amplitudes have the same phase. If the approximant 
is big enough, the inequivalent amplitudes should be nearly 
the same. 

The deviations from symmetry always grow with the 
perp-space reciprocal lattice vector g^. Meanwhile, the in- 
tensities decrease with g^. So the simulation cell should 
preferably be big enough so that, with increasing g-*-, the 
intensities become unmeasurable (or lost in the diffuse scat- 
tering) before the deviations from icosahedral symmetry be- 
come overwhelming. But often the simulation cell must be 
small, for the technical reasons we turn to next. 

5.2.2 Structure factor sum 

In all cases, one must re-evaluate structure factors repeatedly 
— in the case of mode (iii), at every step of the iteration. Even 
if the cell has relatively few tiles, it still has more atoms than 
any unit cell of a metal crystal. (For example, for the canon- 
ical cell tiling, the smallest reasonable approximant is the 
3/2 cubic cell, which contains typically 32 tile nodes (112 
tiles), but some 2500 atoms.) Summing the exponentials is 
relatively expensive in computer time, and the Fast Fourier 
Transform appears to be useless since the sites do not lie on 
a simple grid. Thus, the simulation cell size may be limited 
by the diffraction sum. 

One can save work by partially factorizing the sums. 
Let the index /i label the type of tile (or tile-object), and 
let the (many) allowed orientations of each tile-object be in- 
dexed by Lo. Let be the number of site type on tile-object 



/i, according to the decoration rule in use, and let v be the 
index labeling them; also let M^^ be the multiplicity of sites 
of type ^.v on tile-object fi. Finally, let N^^ be the number 
of tile-objects of type /i and orientation cj in the tiling. Then 
we can define the "tile structure factors" as 

where R^t^j is the reference point of the j-th tile-object of 
type /i in orientation lo. 

Also, we can define the "decoration structure factors" 

as 

F'^'"'°^,y^ = fa{f,u){(l) exp[iq • u^^^,] (27) 

1=1 

Here Uf^^i is the displacement of the decorating atom from 
the reference point of that tile-object, a{iiv) is the atom 
species in orbit jiv, and fa (q) is the form factor for species 
a. Note that F'^^'^° for different orientations to is trivially 
obtained from some reference orientation simply by apply- 
ing the corresponding rotation or reflection to q. This sym- 
metry can reduce the labor, though, only if that rotation or 
reflection takes one of the allowed q vectors to another one. 
The total structure factor then takes the form 

F(q) = F'''%M)F'''''\uM)- (28) 

In a mode (i) process, the first factor remains fixed 
in every term; in mode (ii), the second factor remains fixed. 
Furthermore, it is easy to compute the change in {F^^^'^^^} 
due to each Monte Carlo update move, and this usually af- 
fects relatively few of the different ^uj types. Similarly, if 
we change a parameter in the decoration rule, we only need 
to update F'^°'^° for the orbit ^lv which that parameter 
applies to. 



6 Decagonals 

Decagonal materials present special problems. The 
greatest problem is the seductiveness of apparent two- 
dimensionality, which makes it so much easier to visualize 
the tile packings than in icosahedral cases. In fact, however, 
the stacking is crucial. Physically, it must amplify interac- 
tions within the layers and thus make the structure more 
ordered than a two-dimensional model could be (and per- 
haps more ordered than analogous icosahedral models are). 
In imaging of the structure, averaging may produce highly 
symmetrical patterns even when these are not really present 
in individual layers. 

6.1 Stacking randomness 

First of all, the equilibrium random tiling phase of a 
decagonal cannot be modeled by a random two-dimensional 



1 The same considerations would apply to the other stacked qua- 
sicrystals (octagonal and dodecagonal). 
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tiling The entropy, which is supposed to stabihze 

such a phase, would be proportional to the two-dimensional 
extent of the system and (in the thermodynamic limit) would 
be negligible compared to the volume. (A corollary of this is 
that simulations of decagonals with simulation cells just one 
lattice constant thick in the c direction do not give a valid 
picture of the tiling randomness. However they are valid for 
discovering which atom arrangements are energetically fa- 
vored, and which tiling is appropriate to model that compo- 
sition.) 

Instead, one must permit stacking randomness |^ that 
is, the tiling configuration is similar but not identical from 
one layer to the next. The difference between adjacent lay- 
ers is constrained by rules exactly analogous to the pack- 
ing rules that determine how tiles may adjoin within a layer. 
Unfortunately, stacking randomness seems to be the least 
understood aspect of decagonal structures, both experimen- 
tally [ |60| , |6l| ] and theoretically. [ ]6^ Like everything else in 
decagonals, the stacking rules might depend sensitively on 
composition. 



6.2 Images 

We repeat an old warning about high-resolution transmis- 
sion electron microscope (HRTEM) images. When a sample 
is thick enough that a typical vertical section crosses a stack- 
ing change, then the image becomes (roughly speaking, and 
ignoring the dynamical effects) a two-dimensional projec- 
tion of the scattering density. That projection is always per- 
fectly quasiperiodic, and can not reveal the randomness of 
individual layers. [[l|, pTj ]. 

For this reason, we implore the practitioners of 
HRTEM (and other techniques that project along the di- 
rection of the electron beam, such as the high-angle annu- 
lar dark-field method discussed elsewhere in this volume): 
please try to determine the thickness of the crystal! The per- 
fection of the image can be evidence for the perfection of the 
quasicrystal only to the extent that an upper bound is placed 
on the thickness. 

Multi-layer simulations of a realistic model of 
i(AlNiCo) are feasible [ p^ but not yet realized. We pro- 
pose that a superposition of time steps, as shown in Fig. || 
is a good ersatz for a superposition of layers. The reason is 
that the 2D tilings in adjacent layers should differ by ele- 
mentary tile flips in isolated places, and the same thing us 
true for configurations of time-evolving 2D tilings separated 
by small steps. (The main qualitative difference is that the 
superposition of layers is less random, at long wavelengths, 
than that of time steps, as discussed in Sec. 6.6 of Ref. [jl|].) 

The point that Fig. [s] is intended to illustrate is that 
clusters of (near) 10-fold symmetry emerge from the pro- 
jection which are not present in any layer. Furthermore, 
to accomplish this, one does not need drastic differences 



^2 We used the term "randomness" rather than "disorder", which 
would connote deviation from a particular, ideal structure. In the ran- 
dom tiling case, the ideal structure is inherently random. The stack- 
ing randomness is that a layer disagrees with the adjacent layers, not 
with an ideal layer. We do not have in mind stacking disorder of the 
Hendricks-Teller kind, e.g. adjacent layers of type A when the normal 
pattern alternates type A and type B layers. 



from one layer to the next. This is intended to be a caution 
about structure models derived from electron microscopy, 
by assuming that the image can be interpreted as a two- 
dimensional structure (made of two layers). Even more so, 
it is a caution that the symmetrical ring motifs, which are 
so striking visually in the images, do not necessarily corre- 
spond to an atomic motif. 



6.3 Supertilings 

A final caution about interpreting HRTEM patterns will be 
obvious to most of the microscopists, but perhaps not to oth- 
ers. Supertiles (see Sec. Q) are much commoner in decagonal 
than in icosahedral phases. |^ Quite typically - in analyzing 
experimental images or simulated configurations [ p^ ^3| ] - 
the structure is manifestly a packing of certain small tiles 
(or clusters). Yet they are highly constrained, so one sus- 
pects that a description in terms of a supertiling would be 
more economical. But which one is correct? From initial 
data, perhaps, one can conjecture that certain configurations 
are forbidden, and we can show this implies that the allowed 
configurations are exactly decorations of a certain simple 
supertiling. Yet we then notice (in larger simulation cells) 
one or two exceptions, indicating that the real rule is more 
complicated. (Although wrong, the simple supertiling is not 
useless - it might be quite adequate for a fit of the available 
diffraction data.) 

This observation is intended as a second caution about 
reducing HRTEM images to tilings by placing nodes at the 
centers of pentagonal or decagonal rings of spots. In the ac- 
tual images, there is often a continuum of patterns between 
(say) symmetrical rings and distorted rings: one must almost 
arbitrarily include some nodes and exclude others. Small 
differences in this rule, or errors in its application, may in- 
duce major changes in the rules of the inferred tiling, in the 
sense that a previously forbidden local pattern may start to 
appear, or vice versa. 



7 Conclusion 

All steps of the procedure for fitting a random-tiling struc- 
ture are now available in principle, but they have not yet 
been put into practice. Even though a new kind of "direct 
method" may determine the phases of Bragg amplitudes 
(Sec. this merely produces the ensemble-averaged scat- 
tering density; there remains the harder task of modeling the 
random-tiling ensemble which has been averaged. In accom- 
plishing this, it seems hard to disentangle energy modeling 
from structure modeling. Indeed we suspect that the best 
structure fits of the future will combine total energy calcula- 
tions and diffraction refinements in some fashion. 

In principle, a simple perp-space Debye -Waller fac- 
tor should not suffice to model the relation between the 
ideal and random tilings, but in practice this seemed to work 



^3 One reason may be that the natural inflation scale for a decagonal 



only r = (VS -|- l)/2 cnmp; 



which as mentioned in Sec. 3.1.1 



red to T for a simple icosahedral, 
is often appropriate in models of 



face-centered icosahedral structures. 
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rather well for the rhombohedron tiling (Sec. ||). More study 
is needed: it would be interesting, for example, to apply the 
"minimum-charge" method of Sec. ^ to the Monte Carlo 
"data" of Sec. |[ And if one is given only the random-tiling 
data, what sort of ideal tiling will be reconstructed from it if 
we adopt the factorization approximation and divide out the 
best-fit Debye- Waller factors? 

Decagonal structures are deceptive. They are easy to 
visualize using two-dimensional images, but the equilibrium 
random ensemble is never represented by a two-dimensional 
tiling. It should be a priority - both in experiment and in 
total-energy modeling - to understand the stacking random- 
ness in decagonals. 

At the very start of this paper, we demanded a model 
which contains information beyond the averaged density: 
implicitly that means information about correlations. Yet we 
set "rules of the game" which restricted us to using only 
Bragg data, discarding the diffuse intensity which is in fact 
the Fourier transform of the correlation function. So, an in- 
teresting question for the future is whether there is some sys- 
tematic way to incorporate diffuse data from all of reciprocal 
space (not merely wings of Bragg peaks) in a quantitative fit 
of parameters of a random tiling structure model. 
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Fig. 1. Average electron density p(r) of i(AlPdMn) in physical space, reconstructed from 503 experimental x-ray structure 
factors of ||l6| using signs obtained by the minimum charge algorithm of Sec. g. The plane of the plot is an icosahedral mirror 
plane (13A x ISA" ) centered on a pseudo-Mackay cluster (no atomic surface). Charge density is represented by darkness in 
(a) and by height in (b). There are several instances of atoms (or partial atoms) at unphysically short separation. 




Fig. 2. The same electron density plotted in Figure [l| but for the 5-fold periodic plane. The figure box is a unit cells of this 
plane, so the physical and perp space directions run obliquely The no surface is centered at (0, 0) and (1/2, 1/2), ni at 
(1/2, 0) and (0, 1/2), and bci at (1/4, 1/4) and (3/4, 3/4). A inflation has been applied to both plots to make the 
rounding of features in parallel and perp-space comparable in magnitude (i.e. so the respective Debye-Waller factors are 
similar); this is seen most clearly in the contour plot (a). The height plot (b) shows that the density profile of the no "surface" 
even in the perp-direction, is rounded to an extent comparable to its entire width. 



16 Structure determinations for random-tiling quasicrystals 



Fig. 3. Cross section p(h ) of atomic surface, for ideal and random 
rhombohedron tilings. The horizontal axis is a one-dimensional section 
along the twofold direction; the vertical scale is normalized so the in- 
tegral of the density J d'^h'^ p[h'^), approximated as spherically sym- 
metric, is the same. The rectangular and the flat-topped curves are the 
ideal and random rhombohedron tilings, respectively. The very broad 
Gaussian is a scan, through the bci atomic surface reconstructed from 
Boudard's data | |lq | using the minimum-charge method of Sec. Q (The 
horizontal scale is converted to units of rhombohedron edges assumina 
the model of Ref. |38r.) This last curve differs from a section of Fig. E 
in that only 300 reflections are used here. 
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Fig. 4. Separate Bragg and diffuse diffraction along tlie 2-fold axis of 
reciprocal space for the rhombohedron random tiling (8/5 approximant). 
Vertical bars are a guide for the eye and indicate the location of the 
Bragg peaks. All reciprocal lattice vectors are included with intensity 
over 10^^. Some peaks are labelled using I 
same figure was presented by Tang in Ref. [ 
the Bragg and diffuse parts. 
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Fig. 5. Diffraction intensities from simulated ideal and random tilings, as a function of perp wavevector. 
Filled circles: Bragg peak intensities /s, ideal from ideal 3D Penrose tiling (3DPT) of rhombohedra. Empty 
circles: Bragg intensities /B.rand from the random 3DPT. Crosses: total intensities Jj, rand from the random 
3DPT. integrate d oy er a sphere centered on the Bragg peak, to be compared with the solid line representing 



eq. 
10 



see Sec. 



3.3. Plot (a) shows diffraction from the tiling nodes only, fitted by c = 0.51, D — 3.0 x 
Plot (b) is from a decoration in which node points have form factor + 1, while the midedge and "body- 
center" points (on axis of each prolate rhombohedron) have form factor —1; this is fitted by c = 0.47 and 



D = 4.7 X 10" 
axis label.) 



in eq. d23b. Note: in this and the next two figures, |g | is called "Q_perp" on the figure 
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Structure determinations for random-tiling quasicrystals 
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Fig. 6. Perp Debye-Waller factor in random tiling. At top |F(g^)| is plotted for the ideal and the random 
tiling, showing the intensity reduction. At bottom, the effective DW factor, 2lVcff from eq. (p^, is plotted 
against \ g'^\'^ (labeled "Q_perp^2" in the fipureV its dependence is linear out to the second node of _F(g^), 
showing the validity of the simple formula (Mq) out to surprisingly large wavevectors. Plot (a) shows diffrac- 
tion from the tiling nodes only; plot (b) is from a decoration in which node points have form factor +1, while 
the midedge and "body-center" points (on the axis of each prolate rhombohedron) have form factor —1. 




10"^ 10"' 10"' lO" 
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Fig. 7. The same data as in Fig. ^a) are plotted here (circles') as a com- 
parison of the integrated intensity with the fit from eg. (M), using the co- 
efficients c and D mentioned in the caption of Fig. pl(a). The crosses rep- 
resent a similar comparison for the decorated model shown in Figs, pi b) 
and |b). 



Fig. 8. Top, a snapshot from one time step in a simulation of the (maximally) random Hexagon-Boat-Star 
tiling. The tiling is decorated with atoms (Ni and Co shown filled), as in the d(AlNiCo) model of Ref. |W9p. 
The tile edges are drawn in only in the left portion. Bottom, a superposition of nine successive time steps; see 
how clusters of tenfold symmetry emerge that are not present in the snapshot. This is a poor man's imitation 
of a superposition of different layers, as viewed in an HRTEM image. 



